% Simulation: portfolio selection

clear all
warning off
clc

tic
load 'emp3.mat'
rep = 1000;
TUN = 0:0.1:1;                                                                                                                                 
L = [60 120];       % training sizes
c = [1 2];          % contraint value for GEC estimator by Fan et al. (2012)    
for s = 1:length(DATA)
    % 0. ready data
    DATAt = DATA{s};
    [T,N] = size(DATAt);
    % 0.1 estimate factor parameters
    muF = mean(FACTOR)';
    covF = cov(FACTOR);
    % 0.2 estimate factor loading parameters
    Y0 = DATAt - repmat(RF,1,N);
    l = ones(T,1);
    B0 = [l FACTOR]\Y0;
    muB = mean(B0,2);
    covB = cov(B0');
    PRT{s+1} = [muB(2:end) covB(2:end,2:end)];
    % 0.3 estimate noise parameter
    E = Y0 - [l FACTOR]*B0;
    covE = cov(E);

    % 1. start simulation    
    parfor j = 1:rep
        % 1.1 generate factors and loadings
        F = mvnrnd(muF,covF,T);
        B = mvnrnd(muB,covB,N)';
        Err = mvnrnd(zeros(N,1),covE,T);
        % 1.2 generate portfolios    
        DATAj = [l F]*B + Err;

        % 2. conduct estimation
        SAVEt = [];
        for i = 1:length(L)
            WL = L(i);
            SR = [];
            for t = 1:(T-WL)/12   % roll one-year at a time
                int = 12*(t-1)+1:12*(t-1)+WL;
                ine = 12*(t-1)+WL+1:12*t+WL;
                Xt = DATAj(int,:);
                Xe = DATAj(ine,:);
                % 2.1 l2-relax 
                W = zeros(N,3);
                for opt = 1:3
                    W(:,opt) = l2relax_ps(Xt, TUN, opt);
                end
                % 2.2 GEC by Fan et al. (2012)
                W2 = zeros(N,length(c));
                for s = 1:length(c)
                    W2(:,s) = gec_fun(Xt, c(s));
                end
                % 2.3 construct various portfolios
                P0 = mean(Xe,2);            % simple averaging
                P1 = Xe*W;                  % l2relax (opt = 1,2,3)
                P2 = Xe*W2;                 % GEC
                % 2.4 evaluate by Sharpe Ratio
                P = [P0 P2 P1];             
                SR(t,:) = mean(P)./std(P);
            end
            SAVEt(i,:) = mean(SR);
        end    

        % 3. save results
        SAVEj{j} = SAVEt;
    end
    SAVE{s} = SAVEj;    
end
toc
% 3.1 print variable statistics
PRT{1} = [muF covF];
prt_tab(cell2mat(reshape(PRT,2,2)'))


% 4. generate final results
R = [];
for s = 1:length(SAVE)
    SAVEt = SAVE{s};
    B = length(SAVEt);
    for i = 1:B
        TAR1(i,:) = SAVEt{i}(1,:);
        TAR2(i,:) = SAVEt{i}(2,:);
    end
    R = [R;mean(TAR1);mean(TAR2)];
end
VN = {'60','120'};
prt_tab(R,repmat(VN,1,3),3)





